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Abstract — We  propose  a  new  algorithm  for  real-time  esti¬ 
mation  of  instantaneous  heart  rate  (HR)  from  noise-laden 
electrocardiogram  (ECG)  waveforms  typical  of  unstructured, 
ambulatory  field  environments.  The  estimation  of  HR  from 
ECG  waveforms  is  an  indirect  measurement  problem  that 
requires  differencing,  which  invariably  amplifies  high-fre¬ 
quency  noise.  We  circumvented  noise  amplification  by 
considering  the  estimation  of  HR  as  the  solution  of  a 
weighted  regularized  least  squares  problem,  which,  in  addi¬ 
tion,  directly  provided  analytically  based  confidence  intervals 
(CIs)  for  the  estimated  HRs.  To  evaluate  the  performance  of 
the  proposed  algorithm,  we  applied  it  to  simulated  data  and 
to  noise-laden  ECG  records  that  were  collected  during 
helicopter  transport  of  trauma-injured  patients  to  a  trauma 
center.  We  compared  the  proposed  algorithm  with  HR 
estimates  produced  by  a  widely  used  vital-sign  travel  monitor 
and  a  standard  HR  estimation  technique,  followed  by 
postprocessing  with  Kalman  filtering  or  spline  smoothing. 
The  simulation  results  indicated  that  our  algorithm  consis¬ 
tently  produced  more  accurate  HR  estimates,  with  estimation 
errors  as  much  as  67%  smaller  than  those  attained  by  the 
postprocessing  methods,  while  the  results  with  the  field- 
collected  data  showed  that  the  proposed  algorithm  produced 
much  smoother  and  reliable  HR  estimates  than  those 
obtained  by  the  vital-sign  monitor.  Moreover,  the  obtained 
CIs  reflected  the  amount  of  noise  in  the  ECG  recording  and 
could  be  used  to  statistically  quantify  uncertainties  in  the  HR 
estimates.  We  conclude  that  the  proposed  method  is  robust 
to  different  types  of  noise  and  is  particularly  suitable  for  use 
in  ambulatory  environments  where  data  quality  is  notori¬ 
ously  poor. 
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INTRODUCTION 

The  estimation  of  instantaneous  heart  rate  (HR)  is 
one  of  the  most  vexing  problems  in  physiological 
measurements.4,25  The  difficulties  arise  from  the 
exquisite  sensitivity  of  such  estimations  to  the  smallest 
amount  of  noise  present  in  electrocardiogram  (ECG) 
waveform  recordings  from  which  the  HRs  are  esti¬ 
mated.5  Many  sophisticated  methods  have  been  pro¬ 
posed  to  deal  with  this  problem;2,13,15  however,  by  and 
large,  these  are  based  on  empirical  approaches  that  fail 
to  address  the  nature  of  such  sensitivity. 

The  majority  of  HR  estimation  algorithms  consist 
of  two  major  steps:  QRS  detection  and  HR  estimation. 
In  this  two-step  approach,  emphasis  has  been  placed 
primarily  on  the  first  step,  the  detection  of  QRS 
complexes  in  the  ECG  signal.  QRS  detection  generally 
involves  linear  and  nonlinear  transformations  of  the 
raw  ECG  to  enhance  the  QRS  complexes  in  the  signal, 
and  can  be  achieved  by  techniques  that  focus  on  the 
ECG  amplitude,  its  first  and  second  derivatives,  or  by 
using  digital  filters.13  Unfortunately,  with  the  excep¬ 
tion  of  template-matching  techniques,  the  most  simple 
and  easy-to-implement  algorithms  for  detecting  QRS 
complexes  are  inherently  vulnerable  to  uncertainties  in 
the  determination  of  the  exact  location  of  the  peak  in 
the  complex  and  to  noise  whose  frequency  spectra 
overlap  with  those  of  the  complex,  leading  to  the 
under-  or  overcounting  of  QRS  events  and,  hence,  to 
inaccurate  estimations  of  beat  occurrence  times.13 
Thresholding  algorithms  used  to  detect  the  occurrence 
of  heartbeats  can  partially  alleviate  the  miscounting 
problem;  however,  these  algorithms  are  based  on  time 
windows  of  limited  length  and  may  yield  inaccurate 
estimates.31 

After  a  time  series  of  beat  occurrence  times  has  been 
determined,  the  next  step  is  to  estimate  instantaneous 
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HRs  by  computing  the  difference  between  successive 
beat  occurrence  times.  This  step  is  generally  followed 
by  various  postprocessing  methods  to  smooth  and 
reduce  the  resulting  unphysiological  HR  variations. 
One  class  of  such  algorithms  focuses  mainly  on 
removing  large  spikes  in  the  resulting  HRs  caused  by 
either  ectopic  beats  or  missed/false  QRS  detections,3,26 
necessarily  deemphasizing  noise  (or  error)  in  the  exact 
determination  of  the  beat  occurrence  times.  Another 
class  of  algorithms  focuses  on  improving  the  overall 
quality  of  the  HR  estimates  by  combining  postpro¬ 
cessing  methods,  such  as  Kalman  filter,  with  signal 
quality  assessment  techniques  and  additional  physio¬ 
logical  data.11,17  Furthermore,  other  well-established 
smoothing  techniques,  such  as  spline  smoothing,28 
could  also  be  readily  implemented  to  postprocess  noisy 
HRs.  Nevertheless,  to  the  best  of  our  knowledge,  to 
date,  no  attempt  has  been  made  to  elucidate  the 
mathematical  provenance  of  such  unphysiological  HR 
variations  in  the  estimation  of  instantaneous  HR.  In 
this  article,  we  provide  the  mathematical  rationale  for 
the  sensitivity  of  HR  estimates  to  measurement  noise, 
which  is  ubiquitous  in  ambulatory  environments. 

In  partly  unstructured  environments,  such  as  in  the 
case  of  transport  of  trauma  patients  to  a  trauma  cen¬ 
ter,  the  recorded  HRs  are  notoriously  unreliable.  For 
example,  our  automated  post  hoc  analysis  of  vital-sign 
data  of  898  trauma  patients  collected  during  helicopter 
transport  to  a  Level  I  trauma  center  indicates  that  only 
44%  of  the  HRs  are  of  sufficient  quality  to  be  used  for 
automated  decision  support  and  closed-loop  control.31 
And  our  analysis  of  vital-sign  data  collected  from 
soldiers  using  field-wearable  monitors  during  daily 
physical  activity  shows  that  such  physiological  data  are 
also  notoriously  unreliable.18  The  percentage  of  reli¬ 
able  vital-sign  data  is  expected  to  be  even  lower  in 
more  austere,  uncontrolled  environments,  such  as  a 
battlefield  or  a  mass-casualty  event  setting,  where  the 
challenges  of  obtaining  accurate  measurements  are 
exacerbated  by  the  unexpected  and  quickly  changing 
physiological  status  of  the  casualties  and  other  con¬ 
founding  factors.  To  address  these  challenges,  our 
group  developed  physiological  data  qualification 
algorithms  that  automatically  assess  the  reliability  of 
the  major  vital  signs.7,18,31  While  these  algorithms  have 
been  shown  to  match  the  assessments  made  by  human 
experts  and  significantly  improve  the  accuracy  of 
automated  decision-making  algorithms,8,9,23  they  have 
some  shortcomings:  they  are  not  designed  for  real-time 
analysis,  require  the  availability  of  redundant  sensor 
measurements,  and  operate  based  on  a  fixed  set  of 
physiological  types  of  measurements.  Furthermore,  to 
date,  these  algorithms  have  only  been  used  to  label 
physiological  data  as  reliable  or  unreliable,  without 


replacing  the  unreliable  data  with  the  improved  esti¬ 
mates. 

In  this  article,  we  present  a  new  algorithm  that 
resolves  these  shortcomings  and  provides  a  rigorous 
analysis  of  the  HR  estimation  problem.  We  note  that 
HRs  are  indirectly  estimated  and  that  their  estimation 
can  be  cast  as  a  solution  to  a  least  squares  problem.  We 
argue  that  errors  in  HR  estimation  are  caused  by  the 
sensitivity  of  the  least  squares  solution  (or  naive  dif¬ 
ferencing)  to  noise  in  the  data,  and  further  show  that 
these  errors  can  be  eliminated  by  constraining  the 
solution  of  the  least  squares  problem  using  a  well- 
known  regularization  technique  with  a  weighting 
scheme,27  i.e.,  weighted  regularized  least  squares 
(WRLS),  leading  to  more  consistent  and  reliable  HR 
estimates.  The  representation  of  HR  estimation  as  a 
solution  to  a  least  squares  problem  also  allows  us  to 
analytically  compute  statistical  confidence  intervals 
(CIs)  on  the  estimated  rates.  We  present  results  of 
simulated  and  field-collected  data  that  demonstrate  the 
performance  of  our  algorithm  on  very  noisy  signals 
and  compare  its  performance  against  well-established 
techniques,  such  as  Kalman  filtering  and  spline 
smoothing.  We  show  that  the  proposed  algorithm 
consistently  provides  superior  performance  to  com¬ 
monly  observed  types  of  noise  in  ECG  records. 

METHODS 

QRS  Detection  Algorithm 

This  section  describes  the  first  step  of  the  HR  esti¬ 
mation  algorithm:  QRS  detection.  As  mentioned  ear¬ 
lier,  this  component  of  HR  estimation  algorithms  has 
been  well  developed  and  studied.  To  process  raw  ECG 
waveforms  for  QRS  complex  identification,  we  modi¬ 
fied  and  customized  a  well-known  algorithm20  to  fit  the 
needs  of  our  perspective  applications.  Figure  1  shows  a 
flow  diagram  of  the  algorithm  along  with  the  corre¬ 
sponding  inputs  and  outputs  of  the  four  processing 
stages.  In  the  first  stage,  we  used  a  5-  to  25-Hz  But- 
terworth  band-pass  filter  to  eliminate  non-QRS-related 
frequencies,  and  in  the  second  stage,  we  computed  the 
difference  between  each  two  consecutive  points  of  the 
entire  data  stream  to  amplify  the  sharp  slopes  of 
the  QRS  complex.  After  differencing,  we  squared  the 
resulting  signal  to  make  the  ECG  samples  positive  and 
to  amplify  the  high-frequency  components.  Finally,  in 
the  last  stage,  we  used  a  low-pass  filter  to  enhance 
the  fiducial  marks  of  the  QRS  complex  and  imple¬ 
mented  a  self-adaptive  thresholding  method  to  detect 
QRS  peaks,  reject  noise,  discriminate  T-waves,  and 
search  back  for  missed  QRS  complexes  if  a  detection 
was  not  made  within  a  certain  time  interval.20 
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The  top  panel  in  Fig.  1  shows  an  example  of  a  raw 
0.5-min-long  ECG  segment  from  our  trauma  patient 
database  described  later  in  “Study  Data”  section.  As 
can  be  seen,  the  raw  ECG  waveform  was  contaminated 
with  low-frequency  interference  between  0.1  and 
0.3  min,  which  could  be  attributed  to  a  baseline  drift 
due  to  respiration  or  loose  electrode  contact.  The 
identification  of  QRS  complexes  in  such  a  signal  could 
be  problematic  for  trained  human  experts  and  even 
more  challenging  for  automated  algorithms.  Never¬ 
theless,  the  bottom  panel  in  Fig.  1  shows  that  the 
algorithm  produced  clearly  identifiable  peaks,  denoted 
by  solid  circles,  which  were  detected  using  our 
aforementioned  self-adaptive  thresholding  method. 
Accordingly,  the  final  output  of  the  QRS  detection 
algorithm  of  an  ECG  of  length  T  is  a  time  series  of 
N  monotonically  increasing  cumulative  beat  occur¬ 
rence  times  (Rj),  0<R\<R2<  •  •  •  <Ri<  •  •  •  <  Rn  <  T. 
(Notice  that  we  selected  this  QRS  detection  algorithm 
due  to  its  popularity  and  ease  of  use;  however,  any 
other  QRS  detection  algorithm  could  be  used  in 


conjunction  with  our  HR  estimation  algorithm 
described  next.) 

HR  Estimation  Algorithm 

In  contrast  to  the  previous  section,  this  section 
describes  our  unique  contribution  to  the  HR  estima¬ 
tion  problem.  It  addresses  the  second  step  of  the  HR 
estimation  algorithm,  where  we  estimate  instantaneous 
HRs  from  a  time  series  of  cumulative  beat  occurrence 
times  obtained  by  a  QRS  detection  algorithm.  Two 
important  observations  led  to  the  conceptualization 
and  development  of  this  novel  approach  to  HR  esti¬ 
mation.  The  first  was  that  a  time  series  of  monotoni¬ 
cally  increasing  cumulative  beat  occurrence  times 
needs  to  be  differenced  to  obtain  a  sequence  of 
R-R  intervals  (RRIs).  Specifically,  the  sequence  of 
RRIs  is  obtained  by  subtracting  two  successive 
cumulative  beat  occurrence  times,  Rj  —  Rt- 
i.e.,  RRI  =  {Rz  —  ^i,  R3  —  R2-1  •  •  • ,  Rn  —  Rn—  i}-  Notice 
that,  by  definition,  the  RRI  time  series  represents 
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FIGURE  1.  The  four  stages  of  the  QRS  detection  algorithm  along  with  corresponding  sample  inputs  and  outputs  of  each  of  the 
stages.  The  ECG  fiducial  points,  or  beat  occurrence  times,  are  marked  with  solid  circles  in  the  bottom  panel. 
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the  first-order  difference  of  the  cumulative  beat 
occurrence  times, 

RRI  =  A  R  (1) 

where  A  denotes  the  first-order  difference  operator. 
Equation  (1)  reveals  that  the  desired  RRIs  [or  equiv¬ 
alently  HRs,  where  HR  =  60/RRI  in  beats/min  (bpm)] 
are  not  a  directly  measured  quantity;  rather,  they  are 
indirectly  obtained  by  computing  the  difference  of  the 
directly  observable  cumulative  beat  occurrence  times 
R.  This  observation  that  HR  estimation  is  a  differ¬ 
encing  estimation  problem  has  not  been  discussed  in 
the  literature  so  far  and,  per  se ,  represents  an  impor¬ 
tant  insight  into  the  HR  estimation  problem. 

The  second  observation  was  that  computing  the 
difference  of  an  observable  quantity  containing  some 
measurement  noise  (or  errors)  amplifies  the  noise, 
significantly  contaminating  the  computed  differences. 
This  is  because  subtraction,  or  differencing,  acts  as  a 
high-pass  filter  that  amplifies  high-frequency  noise.12 
In  our  case,  the  right-hand  side  of  Eq.  (1)  contains 
measurement  noise  because  cumulative  beat  occur¬ 
rence  times  cannot  be  determined  with  absolute  cer¬ 
tainty  regardless  of  which  QRS  detection  algorithm  is 
used.  In  fact,  the  American  National  Standards  Insti¬ 
tute  (ANSI)  recommends  a  time  window  of  0.15  s  for 
considering  the  synchronization  between  a  QRS 
detection  algorithm  and  a  reference  annotation.1  In 
addition,  and  more  importantly,  any  QRS  detection 
algorithm  eventually  misses  or  overcounts  a  beat, 
adding  more  noise  to  the  determination  of  the  cumu¬ 
lative  beat  occurrence  times.  The  sensitivity  of  such 
determinations  to  the  slightest  misidentification  of 
peaks  in  ECG  waveforms  has  been  well  documented,5 
and  we  argue  that  this  sensitivity  is  significantly 
amplified  in  the  computation  of  differences  A R,  which 
are  inherent  in  the  calculation  of  RRI.  Therefore,  we 
propose  a  different  approach  to  estimate  HR,  which 
avoids  noise  amplification  during  differencing,  thus 
removing  the  primary  source  of  noise  in  RRI  estima¬ 
tion.  While  all  currently  available  techniques4  attempt 
to  deal  with  the  effect  of  noise  amplification,  the  pro¬ 
posed  approach  removes  the  very  cause  of  such 
amplification. 

In  our  approach,  instead  of  directly  taking  the  dif¬ 
ference  A R  of  the  cumulative  beat  occurrence  times  as 
in  Eq.  (1),  we  reformulated  the  problem  with  R  being 
represented  as  the  integration  of  RRI 

R  =  A  •  RRI  +  e,  (2) 


and  estimated  RRI  as  a  solution  to  an  ordinary  least 
squares  (OLS)  problem24 


RRIqls 


{A1 -A)  l-AT 


■R, 


(3) 


where  R  is  an  N  x  1  vector  of  measured  cumulative 
beat  occurrence  times,  RRI  (RRIols)  is  a  Ax  1  vector 
of  the  corresponding  “true”  (OLS  solution)  RRI  val¬ 
ues,  A  denotes  an  N  x  N  lower  triangular  integration 
matrix  with  all  non-zero  elements  equal  to  one,  s  rep¬ 
resents  a  A  x  1  vector  of  measurement  noise  in  R ,  and 
N  is  the  total  number  of  cumulative  beat  occurrence 
times.  Note  that  Eq.  (3)  represents  the  solution  of  a 
general  OLS  estimation6  (i.e.,  when  A  is  not  invertible) 
and  the  right-hand  side  is  equivalent  to  the  difference 
operation  A R  in  Eq.  (1)  when  A~l  is  analytically 
obtainable,  as  in  the  studied  case.  Using  Eq.  (3),  with 
RRI  expressed  in  seconds,  we  computed  HRqls  = 
60/RRIqls  in  bpm. 

Because  the  OLS  solution  is  equivalent  to  naive 
differencing,  the  resulting  RRIqls  still  contained 
amplified  measurement  noise.  Thus,  to  obtain  a 
smoothed  estimation  of  RRI  without  noise  amplifica¬ 
tion,  we  used  the  WRLS  method,22  which  augments 
the  least  squares  cost  function  \\R  —  A  •  RRI ||2  with  a 
weighting  matrix  W  and  a  penalty  term  to  constrain 
the  variability  of  the  solution.  Thus,  the  goal  is  to 
obtain  RRI  that  minimizes 

(p  •  W  •  RRIols  -AW- RRI||2+A2  •  \\L  •  RRI||2) 

— >  min,  (4) 

where  W  denotes  a  diagonal  N  x  N  weighting  matrix, 
where  the  elements  are  either  zeros  (represented  by 
10-5)  for  spike-like  outliers  detected  in  RRIOLS  (and 
HRols)  via  an  impulse  rejection  filter26  or  ones  for 
non-outliers,  L  denotes  a  smoothing  matrix  that  con¬ 
strains  high-frequency  noise  amplification  in  the  RRI 
estimates  and  produces  a  smooth  and  consistent  solu¬ 
tion,  and  X  represents  a  positive  regularization 
parameter,  which  controls  the  tradeoff  between  the  fit 
to  the  data  and  the  smoothness  of  the  solution.  A 
standard  choice  for  L  (and  the  one  used  here)  is  to  use 
a  (N  —  2)  x  N  matrix  representing  a  second-order 
derivative.22  Conversely,  the  regularization  parameter 
X  is  dependent  on  the  signal-to-noise  ratio  in  the  data 
and  can  be  selected  by  numerous  methods,  such  as 
generalized  cross-validation14  or  the  discrepancy  prin¬ 
ciple.19  Because  the  signal-to-noise  ratio  may  vary  in 
different  data  sets,  and  to  improve  generalization,  we 
customized  X  for  every  patient.  Specifically,  starting 
with  X  =  0  (i.e.,  no  regularization),  we  incrementally 
increased  it  until  the  absolute  time  rate  of  change  of 
the  estimated  HRs  dropped  below  a  specified  threshold 
of  4.0  bpm/s,  which  represents  the  average  absolute 
time  rate  of  change  of  HRs  estimated  from  clean  ECG 
segments  in  our  trauma  patients  database.8 

The  minimization  of  Eq.  (4),  representing  the  WRLS 
solution  for  RRI,  can  be  analytically  obtained 
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RRIwrls  =  (WT  -A1  A  -  W+  l2  ■  LJ  ■  L) 

■  WT  ■  AJ  ■  A  ■  W ■  RRIols-  (5) 

Accordingly,  we  computed  the  WRLS  solution  for 
HR  as  HRWRLs  =  60/RRIwrls- 

Confidence  Interval  Calculation 

Finally,  we  computed  the  Cl  for  the  estimated 
HRwrls  through  a  standard  formulation:10 


Cl  =  HRwrls  ±  C/2  *  \/Var(HRwRLs),  (6) 

where  C/2  denotes  a  percentile  of  a  Student’s  t  distri¬ 
bution  with  a  significance  level  of  a  and  Var(HRWRLS) 
represents  the  variance  of  HRwrls.  The  derivation  of 
Var(HRWRLS)  is  described  in  the  Appendix. 

Study  Data 

Simulated  Data 

To  validate  the  WRLS  algorithm,  we  first  tested  and 
compared  its  performance  using  a  pristine,  i.e.,  a  clean, 
noise-free,  2.5-min-long  ECG  segment  selected  from 
our  trauma  database  (see  “Field-Collected  Data”) 
from  which  we  identified  the  cumulative  beat  occur¬ 
rence  times  R  and  estimated  (through  naive  differ¬ 
encing)  the  “ground  truth”  instantaneous  HRs. 
Subsequently,  we  simulated  two  types  of  noise  in  R: 
small  time  shifts  and  missed/added  beats.21  We  simu¬ 
lated  small  time  shifts,  representing  small  errors  in  the 
detection  of  each  peak  in  a  QRS  complex,  by  adding 
randomly  selected  noise  sampled  from  a  uniform  dis¬ 
tribution  in  the  range  [—b  b],  with  b  set  to  0,  0.05,  0.10, 
or  0.15  s.  The  largest  b  of  0.15  s  was  selected  because 
the  ANSI  stipulates  that  a  timing  mismatch  of  150  ms 
or  less  is  considered  insignificant.1  For  the  second  type 
of  noise,  we  simulated  both  missed  and  added  beats  at 
random  points  in  the  beat  occurrence  time  data 
record.5  In  this  simulation,  we  substituted  0,  10,  20,  or 
30%  of  the  original  data  record,  where  each  time  we 
randomly  decided  to  substitute  a  record  by  deleting  or 
adding  a  beat.  In  total,  we  simulated  16  combinations 
of  noise  conditions  (including  one  noise-free  condi¬ 
tion).  In  addition,  to  separately  investigate  the  per¬ 
formance  of  the  algorithm  on  clusters  of  added  beats 
and  on  clusters  of  missed  beats,  we  simulated  six  more 
noise  conditions,  where  in  each  condition  we  either 
added  beats  or  missed  beats  in  10,  20,  or  30%  of  the 
original  data  record.  Note  that  we  chose  to  introduce 
noise  on  R  instead  of  on  ECGs  because  the  former 
allowed  us  to  directly  test  the  HR  estimation  step,  as 
opposed  to  the  detection  of  QRS  complexes. 


We  simulated  each  of  the  22  noise  conditions  (except 
for  the  noise-free  one)  100  times,  and  computed  the 
corresponding  mean  and  standard  deviation  (SD)  of  the 
root  mean  squared  error  (RMSE)  to  assess  the  ability  of 
four  different  methods  (naive  differencing  (i.e.,  OLS), 
Kalman  filtering,  spline  smoothing,  and  WRLS)  to 
reproduce  the  ground-truth  HRs.  More  specifically,  for 
each  method,  we  first  resampled  the  estimated  and  the 
ground-truth  HRs  to  1Hz  via  linear  interpolation,  and 
then  calculated  the  RMSE  of  the  resampled  rates  by 
taking  the  square  root  of  the  mean  squared  differences  of 
these  two  HRs.  We  applied  Wilcoxon  signed-rank  tests 
to  determine  whether  the  RMSEs  were  statistically  sig¬ 
nificantly  different  from  those  obtained  with  the  WRLS 
algorithm.  The  Wilcoxon  signed-rank  test  is  a  non- 
parametric  statistical  hypothesis  test  for  evaluating  two 
related  samples  or  repeated  measurements  on  a  single 
sample.30  It  can  be  used  as  an  alternative  to  the  paired 
Student’s  Ftest  when  the  population  cannot  be  assumed 
to  be  normally  distributed. 

Field-Collected  Data 

To  validate  the  WRLS  algorithm  using  field-collected 
data,  we  tested  the  estimated  instantaneous  HRs  using 
ECG  recordings  selected  from  our  database  containing 
898  trauma  patients.8  The  physiological  time-series  data 
in  this  database  were  collected  from  injured  patients 
during  transport  via  helicopter  ambulance  from  the 
scene  of  injury  to  the  Level  I  unit  at  the  Memorial 
Hermann  Hospital  in  Houston,  TX.  The  time-series 
variables  were  measured  by  Propaq  206EL  vital-sign 
monitors  (Welch  Allyn;  Skaneateles  Falls,  NY),  down¬ 
loaded  to  an  attached  personal  digital  assistant,  and 
ultimately  stored  in  our  database.  The  physiological 
data  include  ECG  waveforms,  sampled  at  182  Hz,  the 
corresponding  monitor-computed  HR,  and  other  vital- 
sign  data  described  elsewhere.8  Patient  attribute  data, 
such  as  demographics,  were  also  collected  via  chart 
review.  Data  were  collected  and  analyzed  with  the 
approval  of  the  local  and  the  U.S.  Army’s  human  sub¬ 
jects  Institutional  Review  Board,  Fort  Detrick,  MD. 
The  duration  of  the  ECG  recordings  varied  for  different 
patients,  with  an  average  length  of  about  25  min / 
patient.  The  Propaq  processed  the  ECG  with  a 
0.5-40  Hz  bandpass  filter.  Visual  inspection  of  the  ECG 
data  revealed  that  they  were  contaminated  with  the 
various  types  of  noise  reported  in  the  literature.13 

RESULTS 

Performance  of  the  Algorithm  on  Simulated  Data 

We  estimated  HRs  in  the  simulated  data,  i.e.,  the 
cumulative  beat  occurrence  times  identified  from 
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pristine  ECG  segments  with  superimposed  noise  (see 
“Study  Data”  section),  using  four  methods:  naive  dif¬ 
ferencing  (i.e.,  OLS),  Kalman  filtering,  spline  smooth¬ 
ing,  and  WRLS.  The  Kalman  filtering17  and  the  spline 
smoothing28  were  used  as  postprocessing  methods  for 
HRs  computed  through  naive  differencing.  For  a  fair 
comparison,  and  to  eliminate  outliers,  we  used  a  similar 
weighting  scheme  in  the  Kalman  filter  as  the  one  used 
in  the  WRLS  algorithm.  Accordingly,  the  Kalman 
filter  measurement  noise  covariance  matrix  M  (which 
was  a  1  x  1  matrix  in  our  case)  was  computed  as 
M  =  Mo  •  exp(l/w2  —  l),  where  M0  was  selected  using 
the  same  criterion  as  the  one  for  X  in  Eq.  (4)  and  w  cor¬ 
responds  to  the  elements  of  the  weighting  matrix  W,  with 
w  =  1  for  non-outliers  and  w  =  10-5  for  outliers.  Fol¬ 
lowing  Li  et  al.17,  we  set  the  Kalman  filter  1  x  1  state 
noise  covariance  matrix  Q  to  0.1.  Similarly,  before  we 
applied  spline  smoothing28  to  the  HRs  computed 
through  naive  differencing,  we  identified  outliers  and 
substituted  them  through  linear  interpolation.  We 
determined  the  smoothing  parameter  using  the  same 
criterion  used  to  select  X. 

Table  1  summarizes  the  results  (mean  (SD),  except 
for  the  noise-free  condition  in  the  first  row)  of  16  com¬ 
binations  of  noise  conditions  and  shows  that  the  WRLS 
algorithm  significantly  depressed  noise  and  consistently 
yielded  smaller  RMSEs  than  each  of  the  other  three 
methods.  The  WRLS  results  were  considerably  better 
than  those  obtained  with  nai've  differencing,  in  particu¬ 
lar  for  larger  noise  levels.  For  example,  for  the  highest 
noise  level  (missed/added  beats  of  30%  and  uniform 
noise  of  [-0.15  0.15]  s),  the  RMSE  for  the  naive  dif¬ 
ferencing  was  81.79  (29.90)  bpm  while  that  for  the 


WRLS  algorithm  was  more  than  one  order  of  magni¬ 
tude  lower  at  8.03  (1.42)  bpm.  With  respect  to  the  two 
postprocessing  methods,  the  WRLS  algorithm  outper¬ 
formed  the  Kalman  filtering  and  the  spline  smoothing  in 
14  out  of  the  16  noise  conditions,  where  in  each  of  these 
14  conditions,  the  ^-values  of  the  Wilcoxon  signed-rank 
tests  were  below  the  0.001  level  of  significance.  For  the 
two  exceptions  where  the  WRLS  algorithm  did  not  have 
the  smallest  RMSE  (the  first  two  rows  in  the  table), 
the  largest  difference  between  the  WRLS’s  error  and 
that  obtained  with  the  spline  smoothing  was  merely 
0.04  bpm. 

Table  2  summarizes  the  simulated  results  [mean 
(SD)],  where  we  separately  added  beats  or  missed  beats 
in  the  data  record,  and  shows  that  the  WRLS  algo¬ 
rithm  significantly  depressed  noise  and  consistently 
yielded  smaller  RMSEs  than  naive  differencing.  With 
respect  to  the  postprocessing  methods,  for  added  beats, 
the  WRLS  algorithm  outperformed  both  the  Kalman 
filtering  and  the  spline  smoothing  methods  for  each  of 
the  three  noise  levels,  with  p-\ alues  for  the  Wilcoxon 
signed-rank  tests  below  the  0.001  level  of  significance. 
As  for  missed  beats,  the  WRLS  algorithm  had  con¬ 
sistently  smaller  RMSEs  than  the  Kalman  filtering 
method  but  marginally  larger  errors  than  the  spline 
smoothing  method,  with  the  largest  difference  of  0.02 
bpm  being  clinically  insignificant. 

Performance  of  the  Algorithm  on  Field-Collected 
ECG  Waveforms 

To  compare  and  contrast  the  performance  of  the 
WRLS  algorithm  on  field-collected  ECG,  we  selected 


TABLE  1.  Comparison  of  the  mean  (SD)  root  mean  squared  error  (RMSE)  of  HRs  estimated  by  the  four  methods  for  the  16 

combinations  of  simulated  noise  conditions. 


Missed/added 

beats 

Uniform  noise 
level  (s) 

Naive  differencing 
RMSE  (bpm) 

Kalman  filtering 
RMSE  (bpm) 

Spline  smoothing 
RMSE  (bpm) 

Weighted  regularized 
least  squares  (WRLS)  RMSE  (bpm) 

0% 

[0.00  0.00] 

N/A 

0.27 

0.19 

0.21 

[-0.05  0.05] 

4.20  (0.26)* 

2.21  (0.12)* 

1.01  (0.1 1)1 

1.05  (0.12) 

[-0.10  0.10] 

8.74  (0.60)* 

3.57  (0.17)* 

2.07  (0.17)* 

1.40  (0.16) 

[-0.15  0.15] 

14.16  (1.15)* 

5.19  (0.33)* 

4.21  (0.43)* 

1.71  (0.21) 

10% 

[0.00  0.00] 

44.72  (19.85)* 

2.48  (0.47)* 

1.21  (0.39)* 

1.04  (0.30) 

[-0.05  0.05] 

46.68  (20.26)* 

3.57  (0.39)* 

2.70  (0.62)* 

2.09  (0.37) 

[-0.10  0.10] 

49.17  (21.06)* 

5.12  (0.67)* 

4.97  (1.03)* 

3.13  (0.52) 

[-0.15  0.15] 

50.21  (21.55)* 

6.67  (0.92)* 

7.36  (1.24)* 

4.66  (0.81) 

20% 

[0.00  0.00] 

69.13  (24.20)* 

3.24  (0.35)* 

2.24  (0.60)* 

1 .85  (0.44) 

[-0.05  0.05] 

67.11  (22.22)* 

4.74  (0.62)* 

4.75  (1.18)* 

3.35  (0.66) 

[-0.10  0.10] 

66.62  (20.90)* 

7.05  (1.19)* 

8.21  (1.62)* 

5.14  (1.01) 

[-0.15  0.15] 

68.19  (22.73)* 

8.47  (1.86)* 

10.39  (2.20)* 

6.67  (1.20) 

30% 

[0.00  0.00] 

80.25  (22.94)* 

3.92  (0.52)* 

3.37  (0.92)* 

2.75  (0.65) 

[-0.05  0.05] 

81.42  (24.96)* 

5.91  (1.10)* 

6.83  (1.78)* 

4.66  (0.98) 

[-0.10  0.10] 

80.81  (22.03)* 

8.34  (1.93)* 

10.68  (2.50)* 

6.90  (1.35) 

[-0.15  0.15] 

81.79  (29.90)* 

10.17  (2.37)* 

13.52  (3.00)* 

8.03  (1.42) 

SD:  standard  deviation;  HR:  heart  rate;  N/A:  not  applicable,  as  naive  differencing  was  used  to  compute  the  ground  truth  HRs  based  on  noise- 
free  condition;  *p-value  <0.001;  tp-value  <0.01. 
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TABLE  2.  Comparison  of  the  mean  (SD)  root  mean  squared  error  (RMSE)  of  HRs  estimated  by  the  four  methods  for  the  six 

simulated  noise  conditions  with  either  added  or  missed  beats. 


Noise  type 

Noise  level 
(%) 

Naive  differencing 
RMSE  (bpm) 

Kalman  filtering 
RMSE  (bpm) 

Spline  smoothing 
RMSE  (bpm) 

Weighted  regularized 
least  squares  (WRLS)  RMSE  (bpm) 

Added  beats 

10 

67.16  (33.12)* 

2.80  (0.21)* 

1.80  (0.38)* 

1.48  (0.26) 

20 

106.05  (31.95)* 

4.21  (0.43)* 

4.05  (0.71)* 

3.06  (0.47) 

30 

134.93  (35.90)* 

8.82  (1.45)* 

10.89  (2.03)* 

7.31  (1.13) 

Missed  beats 

10 

16.18  (0.78)* 

1 .02  (0.33)* 

0.41  (0.06)* 

0.43  (0.07) 

20 

24.52  (0.89)* 

1.76  (0.41)* 

0.61  (0.10)* 

0.63  (0.11) 

30 

32.05  (1.33)* 

2.36  (0.50)* 

0.83  (0.18)* 

0.84  (0.18) 

SD:  standard  deviation;  HR:  heart  rate;  *p-value  <0.001,  Rvalue  >0.05. 
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FIGURE  2.  Example  of  noisy  ECG  signals  collected  from  a  patient  during  transport  to  a  trauma  center  and  the  corresponding 
heart  rates  [HRs;  in  bpm]  estimated  by  four  distinct  methods.  The  top  panel  depicts  an  ECG  segment  contaminated  with  spike 
noise  between  0.4  and  1.2  min  and  electromyographic  noise  between  1. 4-2.5  and  3.7-4.0  min.  The  second  and  third  panels 
compare  HRs  estimated  by  the  proposed  weighted  regularized  least  squares  (WRLS)  algorithm  against  those  obtained  by  Propaq, 
Kalman  filtering,  and  spline  smoothing.  The  bottom  panel  illustrates  the  WRLS-estimated  HRs  along  with  the  corresponding  95% 
confidence  intervals. 


two  different  representative  patients  from  our  trauma 
database  whose  ECGs  were  contaminated  with  com¬ 
monly  observed  noise  artifacts.  Figure  2  shows  4  min 
of  a  normalized  ECG  record  {top  panel  and  inset  dis¬ 
playing  a  smaller  scale)  for  one  patient  and  the  asso¬ 
ciated  HR  estimates  {middle  two  panels)  calculated 


from  a  Propaq  206EL  vital-sign  monitor,  Kalman  fil¬ 
tering,  spline  smoothing,  and  the  proposed  WRLS 
algorithm.  The  ECG  in  Fig.  2  was  contaminated  with 
spike  noise  and  electromyographic  (EMG)  noise.  Such 
spike  noise,  consisting  of  high-frequency  components 
resembling  those  of  QRS  complexes,  is  common  in  our 
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trauma  patient  database,  and  EMG  noise  due  to 
muscle  contractions  is  ubiquitous  in  ECG  signals.13 
Because  EMG  noise  partially  overlaps  with  the  QRS 
spectra  in  the  frequency  domain,  most  HR  monitors, 
including  Propaq,  use  some  kind  of  EMG  noise  fil¬ 
tering.  However,  the  complete  removal  of  this  kind  of 
contaminant  is  usually  not  possible.  In  Fig.  2,  the 
largest  amount  of  EMG  noise  was  observed  between 
1.4-2. 5  and  3. 7-4.0  min,  outside  of  which  the  ECG 
was  relatively  clean,  except  for  the  spike-noise-con¬ 
taminated  segment  between  0.4  and  1.2  min.  The  ECG 
contamination  was  reflected  in  the  Propaq-calculated 
HRs  (, second  panel),  which  showed  very  large  varia¬ 
tions  between  1.4-2. 5  and  3. 7-4.0  min,  changing  from 
100  to  200  bpm  within  a  few  seconds,  and  smaller  but 
marked  variations  between  0.5  and  1.2  min.  Such  rapid 
HR  excursions  are  usually  unphysiological  and  indi¬ 
cate  the  inability  of  the  algorithm  to  provide  continu¬ 
ously  smooth  point  estimates.  It  is  not  clear  why  the 
Propaq-calculated  HRs  also  showed  marked  variations 
between  2.5  and  2.8  min,  when  the  ECG  recording  was 
relatively  clean.  In  contrast,  the  WRLS  algorithm 
provided  smooth  and  consistent  estimates  for  the 


whole  recording,  indicating  its  ability  to  handle  both 
types  of  noise.  Note  that  during  the  periods  of  clean 
ECG  segments,  the  two  algorithms  yielded  relatively 
similar  HR  estimates.  The  third  panel  in  Fig.  2  com¬ 
pares  HR  estimates  obtained  by  the  WRLS  algorithm 
with  those  from  the  Kalman  filtering  and  the  spline 
smoothing  methods.  While  overall  there  was  good 
agreement  among  the  three  HR  estimates,  the  Kalman 
filter  tended  to  fit  a  flat  line  no  matter  whether  the 
“intrinsic”  HR  was  changing  or  not,  and,  in  compar¬ 
ison  with  our  algorithm,  the  spline  smoothing  exhib¬ 
ited  small  overestimations,  for  example,  at  0.5,  0.7,  and 
0.9  min. 

Figure  3  shows  similar  comparisons  for  three  addi¬ 
tional  types  of  commonly  found  noise  in  ECG 
recordings:  motion  artifacts,  electrode  contact  noise, 
and  instrumentation  saturation  noise.  Motion  artifacts 
are  transient  (but  not  step)  baseline  variations  caused 
by  changes  in  the  electrode-skin  impedance  due  to 
electrode  motion13  and  are  commonly  observed  during 
patient  transport.  Electrode  contact  noise  is  generally 
caused  by  a  loose  contact  between  an  electrode  and  a 
patient’s  skin,  and  saturation  noise  is  usually  caused  by 


Electrode  contact  +  Instrumentation  saturation  noise 
Motion  artifact  _ h. 
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FIGURE  3.  Example  of  noisy  ECG  signals  collected  from  a  patient  during  transport  to  a  trauma  center  and  the  corresponding 
heart  rates  [HRs;  in  bpm]  estimated  by  four  distinct  methods.  The  top  panel  depicts  an  ECG  segment  contaminated  with  electrode 
contact  and  instrumentation  saturation  noise  between  1.5  and  2.5  min  and  motion  artifact  between  0.6  and  0.7  min.  The  second  and 
third  panels  compare  HRs  estimated  by  the  proposed  weighted  regularized  least  squares  (WRLS)  algorithm  against  those  obtained 
by  Propaq,  Kalman  filtering,  and  spline  smoothing.  The  bottom  panel  illustrates  the  WRLS-estimated  HRs  along  with  the  corre¬ 
sponding  95%  confidence  intervals. 
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amplification  of  artifacts  exceeding  the  data  range.13 
The  latter  two  types  of  noise  significantly  affected  the 
ECG  between  1.5  and  2.5  min,  causing  Propaq’s  esti¬ 
mates  to  change  from  around  100  bpm  to  as  low  as 
70  bpm  in  just  a  few  seconds.  Between  0.6  and  0.7  min 
of  the  recording,  the  Propaq-estimated  HRs  displayed 
a  “bump”  above  the  100  bpm  rate,  reflecting  the  effects 
of  the  ECG  motion  artifact  in  its  estimates.  In  con¬ 
trast,  the  WRLS  algorithm  produced  smooth  and 
consistent  estimates.  Notice  that  the  ECG  was  rela¬ 
tively  clean  before  0.5  min  and  after  2.5  min  of  the 
recording,  during  which  periods  the  Propaq’s  and  the 
WRLS’s  estimates  practically  coincided  with  each 
other,  indicating  that  during  the  clean  ECG  segments 
both  algorithms  worked  correctly.  When  compared 
with  the  WRLS  algorithm,  the  Kalman  filtering  and 
the  spline  smoothing-estimated  HRs  showed  a  similar 
pattern  ( third  panel),  although  the  differences  were 
more  pronounced  than  those  in  Fig.  2.  In  particular, 
during  the  clean  ECG  recording  before  0.5  min  the 
Kalman  filtered  HRs  displayed  a  ~5%  deviation  from 
those  obtained  with  the  other  two  methods,  and  at  ~1 .7 
min  the  spline  smoothing  overestimated  the  WRLS  by 
~3%  and  the  Kalman  filter  by  ~8%. 

Confidence  Intervals 

The  bottom  panels  in  Figs.  2  and  3  illustrate  the 
WRLS  estimated  HRs  and  their  corresponding  95% 
CIs.  The  width  of  the  CIs  should  reflect  the  uncertainty 
in  the  HR  estimates,  with  smaller  width  indicating 
lower  uncertainty,  and  vice  versa.  Accordingly,  ECG 
segments  containing  larger  amounts  of  noise  artifacts 
should  yield  wider  CIs.  For  example,  the  estimated  CIs 
in  the  bottom  panel  in  Fig.  2  were,  as  expected,  large 
between  0.4-1. 2  and  1.4-2. 5  min,  indicating  uncer¬ 
tainties  in  the  HR  estimates  associated  with  spike  and 
EMG  noises.  Outside  of  these  regions,  the  CIs  were 
narrow,  reflecting  small  uncertainties  in  the  HR  esti¬ 
mates  because  of  the  good  quality  of  the  ECG  signal. 
Figure  3  shows  a  similar  pattern;  the  CIs  were  wider 
for  noisy  regions  in  the  ECG  and  narrower  for  good- 
quality  ones. 


DISCUSSIONS  AND  CONCLUSIONS 

This  article  provides  a  formal  mathematical  for¬ 
mulation  to  compute  instantaneous  HR  in  real  time 
from  noise-laden  ECG  signals  typical  of  field  envi¬ 
ronments.  We  noted  that  HR  estimation  is  an  indirect 
measurement  problem  that  requires  the  computation 
of  differences  between  successive  beat  occurrence 
times,  which  invariably  amplifies  noise  (i.e.,  the 
uncertainties  in  the  determination  of  the  beat 


occurrence  times).  We  circumvented  noise  amplification 
by  casting  the  estimation  of  HRs,  from  beat  occurrence 
times,  as  the  solution  of  a  weighted  regularized  least 
squares  problem,  thereby  eliminating  the  need  for 
post  hoc  smoothing  of  HR  signals  corrupted  with 
amplified  noise.  Hence,  the  proposed  formula¬ 
tion  is  fundamentally  different  from  existing  tech¬ 
niques3,4,11,17,26  in  that  while  these  solutions  attempt  to 
smooth  noise-amplified  HRs  with  poor  signal-to-noise 
ratios,  the  proposed  algorithm  avoids  noise  amplifica¬ 
tion  altogether. 

A  unique  advantage  of  the  proposed  WRLS  algo¬ 
rithm  is  that  the  resulting  solution  is  optimal  in  the  least 
squares  sense  and,  hence,  it  yields  a  closer  estimate  to 
the  true  (unknown)  solution  than  the  standard  least 
squares  approach.16  Therefore,  in  principle,  the  method 
of  regularized  least  squares  can  improve  generalization 
with  respect  to  different  data  sets.  However,  regulari¬ 
zation  is  inherently  biased,  as  it  trades  off  a  solution 
with  a  larger  variance  and  a  smaller  bias  for  one  with  a 
smaller  variance  but  a  slightly  larger  bias,  where  the 
bias  is  controlled  by  the  selection  of  the  smoothing 
matrix  L  and  the  regularization  parameter  X.  Hence,  in 
practice,  regularization  can  also  limit  generalization. 
To  avoid  this  potential  problem,  we  set  L  as  a  second- 
order  derivative,  which  is  the  most  general  choice  for 
imposing  smoothness  in  the  estimates  of  R-R  inter¬ 
vals,22  and  customized  X  for  every  patient  to  allow  the 
estimated  HRs  to  vary  as  much  as  possible  within  the 
“normal”  range  (i.e.,  ±4.0  bpm/s).  We  also  introduced 
a  weighting  scheme  to  eliminate  the  effects  of  outliers, 
which  could  otherwise  significantly  bias  the  results  even 
with  the  use  of  regularization. 

Another  advantage  of  the  proposed  algorithm  is  the 
ability  to  infer  a  measure  of  uncertainty  of  the  esti¬ 
mated  HRs  in  the  form  of  CIs.  This  is  achieved  by 
casting  HR  estimation  as  a  solution  to  a  least  squares 
problem,  directly  yielding  analytic,  statistically  based 
CIs  for  the  estimated  rates  that  account  for  noise  in  the 
ECG,  a  feature  that  is  lacking  in  any  prevailing  algo¬ 
rithm.  This  is  in  contrast  with  our  previous  attempts  to 
qualify  physiological  measurements  based  on  empirical 
rules  with  no  statistical  justification.7,18,31  Moreover, 
our  method  is  generic,  requiring  only  the  availability  of 
a  single  waveform  signal.  This  independence  of  other 
physiological  measurements  promotes  usability  across 
different  waveform  signals,  e.g.,  respiratory  wave¬ 
forms,  and  portability  across  vital-sign-monitoring 
devices  and  decision-support  systems. 

To  validate  our  WRLS  algorithm,  we  assessed  its 
performance  with  both  simulated  and  field-collected 
data  against  established  HR  estimation  methods.  For 
the  simulated  data,  the  WRLS  algorithm  yielded  esti¬ 
mation  errors  (i.e.,  RMSEs)  that  were  more  than  one 
order  of  magnitude  smaller  than  those  obtained 
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through  naive  differencing  (i.e.,  OLS).  When  we  post- 
processed  these  HRs  with  two  different  methods 
(Kalman  filter  and  spline)  to  smooth  the  estimated 
rates,  we  obtained  significantly  improved  results. 
However,  the  attained  RMSEs  were  still  higher  by  as 
much  as  67%  than  those  obtained  with  the  WRLS 
algorithm,  which  consistently  yielded  superior  results. 
This  supports  our  assertion  that  rather  than  attempt¬ 
ing  to  smooth  noise-amplified  HRs  with  poor  signal- 
to-noise  ratios,  one  should  avoid  noise  amplification 
altogether.  When  we  applied  the  WRLS  algorithm  to 
the  898  patients  in  our  trauma  database,  it  reliably 
produced  improved  HR  estimates  that  showed 
observable  differences  when  compared  with  those  pro¬ 
vided  by  the  monitoring  device  (i.e.,  Propaq)  and  the  two 
postprocessing  algorithms.  In  the  presence  of  noise  fre¬ 
quently  observed  in  field-collected  ECG,  the  Propaq 
often  yielded  HRs  that  exhibited  very  rapid  excursions, 
suggesting  that  the  estimated  HRs  were  unphysiologi- 
cal.  In  contrast,  the  Kalman  filter  failed  to  express  any 
variability  in  its  smoothing,  as  it  systematically  fitted  a 
flat  line  regardless  the  amount  of  noise  in  the  naive  dif¬ 
ferencing  HR  estimates.  These  estimates  could  be  im¬ 
proved,  however,  by  jointly  employing  signal  quality 
indices  and  using  additional  vital-sign  information.17 
Conversely,  spline  smoothing  frequently  exhibited  lar¬ 
ger  variations  in  the  HR  estimates,  which  we  attribute  to 
noise  amplification.  In  addition,  the  CIs  estimated  by  the 
WRLS  algorithm  correctly  reflected  the  uncertainties  in 
the  HR  estimates  with  respect  to  noise  contamination  in 
the  ECG  recordings. 

One  potential  limitation  of  the  proposed  method  is 
the  inability  to  track  true,  sudden  large  changes  in  HR. 
Because  we  imposed  a  constraint  of  4.0  bpm/s  on  the 
maximum  time  rate  of  change  of  “normal”  HRs,  the 
proposed  method  could  compromise  the  accuracy  of 
the  instantaneous  HR  estimates  if  the  true  HR  changes 
were  larger  than  4.0  bpm/s.  In  this  case,  the  method 
would  correctly  detect  the  onset  of  elevations  or 
reductions  in  HR  without  time-lags;  however,  it  would 
produce  a  delay  in  reaching  the  actual  elevated  or 
reduced  HR  value.  Such  a  delay  would  be  proportional 
to  the  difference  between  the  true  absolute  time  rate  of 
change  of  the  HR  and  the  imposed  constraint  of 
4.0  bpm/s.  Naturally,  this  limitation  could  be  allevi¬ 
ated  by  arbitrarily  increasing  the  constraint;  however, 
setting  the  constraint  too  high  would  also  diminish  the 
accuracy  of  the  algorithm  by  including  excessive  noise 
in  the  estimated  HRs. 

As  biosensors  become  ubiquitous  in  everyday  life,  it 
is  important  that  we  continue  improving  algorithms 
for  the  real-time  estimation  of  vital  signs.  For  both 
civilian  and  military  applications,  it  is  particularly 
important  to  infer  reliable  values  for  HRs — arguably 


one  of  the  most-used  vital  signs — collected  from  aus¬ 
tere,  unstructured  environments,  such  as  a  battlefield, 
during  the  transport  of  trauma  patients,  in-home  care 
of  elderly  patients,  and  in  the  monitoring  of  active 
individuals  during  physical  activity,  where  the  original 
ECG  data  are  prone  to  be  contaminated  with  noise 
artifacts.  The  study  proposed  here  is  a  step  in  this 
direction,  as  it  allows  for  more  consistent  estimation  of 
HRs  and  of  other  ECG  features,  such  as  HR  vari¬ 
ability,  which  are  highly  sensitive  to  ECG  noise. 


APPENDIX 

Because  HRwrls  is  reciprocal  to  RRIwrls,  we 
approximated  Var(HRWRLS)  using  the  mean  and  var¬ 
iance  of  RRIwrls  through  a  Taylor  expansion,  i.e., 

Yar  (HRwrls) 

*  h - - -A  Var(RRIwRLs),  (Al) 

\Jmean(RRIWRLS)]  J 

where  Var(RRIWRLS)  was  estimated  as  the  diagonal  of 
the  covariance  of  RRIwrls,  Cov(RRIwrls),  which 
was  calculated  from  Eqs.  (2),  (3),  and  (5)  as 

Cov( RRIwrls)  =  C-B-  Co v(e)  •  BT  •  CT,  (A2) 

where  B  =  (AT  ■  A)~l-AT,  C  =  (W1  ■  AJ  ■  A  ■  W  +  A2- 
LJ  ■  L)~l  ■  WT  ■  At  ■  A  ■  W,  and  Cov(e)  denotes  the 
covariance  of  the  measurement  noise,  which  was  esti¬ 
mated  as  a  diagonal  matrix  with  elements  equal  to  half 
of  the  square  of  the  residual  between  RRIqls  and 
RRIwrls-29 
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